Mesoscopic simulations of the counterion-induced electroosmotic flow - a comparative 

study 
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We present mesoscopic simulations ol the counterion-induced electroosmotic flow in different elec- 
trostatic coupling regimes. Two simulation methods are compared, Dissipative Particle Dynamics 
(DPD) and coupled Lattice- Boltzmann/Molecular Dynamics (LB/MD). A general mapping scheme 
to match DPD to LB/MD is developed. For the weak-coupling regime, analytic expressions for 
the flow profiles in the presence of partial-slip as well as no-slip boundary conditions are derived 
from the Poisson-Boltzmann and Stokes equations, which are in good agreement with the numerical 
results. The influence of electrofriction and partial slip on the flow profiles is discussed. 



INTRODUCTION 

Microfluidic devices like bio-MEMS (micro- 
electronical-mechanical-systems) and bio-NEMS (nano- 
electronical-mechanical-systems) have attracted broad 
interest over the last years due to their huge potential in 
biotechnology [l|, [3| • The flow profiles in such micro- or 
nanosized devices are strongly influenced by the proper- 
ties of the boundaries due to the large surface-to- volume 
ratio in these systems. Surface characteristics like 
the wetting behavior and/or slippage have a dramatic 
effect on the microscopic flow, leading to sometimes 
unexpected behavior 0]. 

One particularly important mechanism is electroos- 
motic transport: in contact with a liquid, many materials 
commonly used in nanotechnology {e.g., polydimethyl- 
siloxane (PDMS)) become charged due to ionizations of 
surface groups [J|. As a consequence, surfaces are often 
covered by a compensating counterion layer In an 
external electric field, the ions are driven in one direc- 
tion, dragging the surrounding solvent with them. As a 
result, a flow is induced in the fluid, the electroosmotic 
flow (EOF). This electrokinetic effect has numerous con- 
sequences. For example, it alters drastically the migra- 
tion dynamics of mesoscopic objects like polyelectrolytes 
or colloids Q. In microchannels, the EOF generated 
at the channel walls results in a total net flow, which 
is technologically attractive because it can be controlled 
and manipulated more easily on the submicrometer scale 
than pressure- or shear-driven flow. 

As analytical predictions for flow in such complex sys- 
tems are often hard to derive 0|, numerical simulations 
are used to investigate it in detail. The requirements 
on computer simulations are high: They must include 
the full hydrodynamics of the fluid while guarantee- 
ing, at the same time, optimal computational efficiency. 
This has led to the development of a number of coarse- 
grained mesoscopic simulation schemes. Several methods 
like Lattice Gas Automata the Lattice-Boltzmann 



Method (LB) p, p , UJj^ , Dissipative Particle Dynamics 
(DPD) [lil,|mll|7and Multi-Particle Colhsion Dynam- 
ics (MPC) [IJ, ll5[ are nowadays used to model liquids 
at a coarse-grained level. A completely discretized ap- 
proach to study EOF via a generalized LB method has 
been proposed by Warren [l6|, and later extended by 
Capuani and coworkers [l7j . relying on a solution of the 
electrokinetic equations on the LB lattice. 

Compared to atomistic Molecular Dynamics simu- 
lations, these approaches give access to much longer 
time- and length scales [18[ and are therefore suited 
to study the long-time behavior of soft matter systems 
and transport phenomena. Although the theoretical 
background of these methods is well understood, their 
lattice/ofF-lattice and thermal/athermal character im- 
pedes a straightforward mapping between them. Com- 
parative studies of specific soft matter problems with 
different simulation methods are therefore scarce. In 
this paper we present the results of DPD- and coupled 
Lattice-Boltzmann/Molecular Dynamics (LB/MD) simu- 
lations for the counterion-induced EOF without salt ions 
in a slit pore. Although this geometry does not allow 
to test the full coupling between ions and solvent, it has 
the advantage of possessing an analytical solution, which 
allows a precise test of the accuracy of the methods. A 
systematic approach to match DPD and LB/MD simu- 
lations is presented. 

As mentioned earlier, flow profiles depend heavily on 
the boundary conditions at the surfaces [1]. We de- 
rive an analytic equation for the counterion-induced EOF 
in the presence of no-slip as well as partial-slip bound- 
aries for the regime of "weak electrostatic coupling" , the 
regime where the Poisson-Boltzmann theory is valid, and 
compare it to our simulation results. Furthermore, we 
also study the influence of the discrete character of wall 
charges, compared to perfectly homogeneous walls. 

The paper is organized as follows. In the flrst section 
we derive the analytical solution for the flow proflle of the 
counterion-induced EOF in the presence of arbitrary slip 
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conditions. The two mesoscopic simulation methods are 
introduced and the general mapping scheme is discussed 
in section III. Section IV gives the details of the simula- 
tions and section V focuses on the numerical results. We 
conclude with a brief summary. 



THEORETICAL BACKGROUND 

When brought in contact with a liquid, most materials 
acquire charges cither by the ionization or dissociation of 
surface groups or the adsorption of ions from solution. In 
microchannels, the walls are thus typically charged and 
attract a layer of compensating counterions. In the fol- 
lowing, we consider the simple situation of fluid flow in 
a slit channel with charged walls and dissolved counteri- 
ons. Additional charges {e.g., salt ions) are not present. 
Depending on the relative strength of Coulomb interac- 
tions in the system, one distinguishes between several 
electrostatic coupling regimes, which can be realized by 
changing the surface charge densities, by using counteri- 
ons of different valency, or by varying the temperature. 
One parameter that discriminates between regimes is the 
electrostatic coupling constant 



(1) 



which gives the strength of electrostatic interactions be- 
tween the countcrion and the surface compared to the 
thermal energy. The parameters entering S are the Bjer- 
rum length Xb — /AnerkBT, the unit charge e, the 
counterion valency Z, the thermal energy ksT, the di- 
electric constant er, and the surface charge density a a 
(19I . [iol . 21 1 . Low values of S correspond to the "weak 
coupling limit" , where the counterions at the wall are 
distributed in a broad diffusive layer which is well de- 
scribed by the Poisson-Boltzmann theory [23, 24|. This 
regime is realized at moderate surface charge densities 
and for monovalent ions. High values of S are obtained 
for high surface charge densities, low temperatures, or 
multivalent charges, and lead to the formation of nearly 
flat, hi ghly adsorbed and massively correlated counterion 
layers |22l |. If in addition the two counterion layers are 
well separated, i.e., the slit width d is much larger than 
the lateral distance between ions a [d/a > 1 [20|), one 
enters the "strong coupling limit" , which in terms of the 
electrostatic coupling constant is given for S > 
with the Gouy-Chapman length /i = (27rZAscrA)~^- 

In the following we focus on the analytical derivation 
of the EOF profiles in the weak coupling regime. The 
electrostatic potential "tpiz) and the ion density distribu- 
tion p{z) can be calculated from the Poisson-Boltzmann 
equation [H,!!!, which reads in our case (slit channel, 
no salt, counterions only) 



9z2 



Ze , , Ze _z£^ 
^{z) = -—p{z) = ~—poe ^bt (2) 



where Ze is the charge of the ions, and po the reference 
ion density at -0 = 0. If an additional electric field 
is applied, the counterions move in the direction of the 
field and drag the solvent particles along. A net flow is 
induced, the counterion-induced EOF. The actual fiow 
profile in the limit of low Reynolds numbers can be cal- 
culated with the Stokes equation, which reads in absence 
of pressure gradients 
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(3) 



with the dynamic viscosity 77. We impose the most gen- 
eral hydrodynamic boundary condition, the partial slip 
condition 



dzVx{z)\., 



--±ZE 



Ty-Wx(±2s) 
(IB 



(4) 



at the hydrodynamic wall boundaries ±zb with the slip 
length 8b. Comparing Eq. ([3]) with Eq. ([2]), we obtain 
by straightforward integration of Eq. ([3]) with @j the 
general relation 



Vx{z) 



(V'(zs) - tP{z) - SBi'\zB)). 



(5) 



Specifically, the solution of Eq. ^ in the slit geometry 
is given by ^ 



kBT 
Ze 



log(cOS'^(Kz)) 



(6) 



with the screening constant ~ {Ze)'^pQ/2erkBT, where 
Pa is the counterion density in the middle of the channel. 
The counterion density distribution is then given by 



and the electric field is 



PO 



cos^(kz) ' 



(7) 



E{z) = -— ^ z = ^ tan . 8 

oz Ze 

For the flow profile, we finally get 



Vx{z) 



-Ex{\og{cos^ [kzb)) 



AttXbZt] 

— log(cos^(Kz)) + 2k6b tan(KZB)), (9) 



where we have expressed in terms of the Bjerrum 
length Xb. 

If the surface charges are high, the ions have high va- 
lency, or the temperature is low, charge correlations and 
charge fiuctuations become important and mean-field ap- 
proaches like the Poisson-Boltzmann approach fail. In 
this case, asymptotic analytical expressions for the ion 
distribution in the channel can be obtained with field- 
theoretic methods [1^, 2^, 21 1. If the plate distance d is 
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much larger than the Gouy-Chapman length, the coun- 
terion density distribution between two planar highly 
charged surfaces is given by 



(1 - e-d^/t") 



This result characterizes the "strong coupling limit" , 
corresponding to two independent highly self-correlated 
counterion layers. 



SIMULATION METHODS 
Dissipative Particle Dynamics 

Dissipative Particle Dynamics (DPD) is a coarse- 
grained particle-based simulation method for fluid dy- 
namics. It is Galilean- invariant, generates a well-defined 
canonical ensemble at equilibrium and guarantees the 
conservation of momentum. The system consists of a set 
of particles with continuous positions and velocities 
whose time evolution is described by Newton's equations 
of motion. The forces on one particle are given by 



^DPD 



(11) 



with the conservative force F,^ , a two-particle dissipative 
interaction 



(12) 



with fij = Tij/vij, Tij = Tj — Ti, and a corresponding 



stochastic force 



(13) 



The weight function Lo{r) can be chosen at will - here we 
use the most common form uj{r) = \ — r/vc for r < 7\. 
{oj{r) = otherwise), with the DPD cut-off radius rc- 
The random numbers Qj are symmetric and otherwise 
uncorrelated {Qj — Qi) with mean zero and variance one. 
The symmetry property F^^ = —FJ^, and F^^ = —F-f 
ensures that the total momentum is conserved in the ab- 
sence of spatially varying external (conservative) poten- 
tials. 

The hydrodynamic boundary condition at the walls (|4]) 
is realized with a recently developed method [2^ that al- 
lows to implement arbitrary partial-slip boundary condi- 
tions: We introduce an additional coordinate-dependent 
viscous force that mimicks the wall/fluid friction 



with a dissipative contribution 

Ff = -7L UJl{z) (Vj - \wall) 



(14) 



(15) 



coupling to the relative velocity (v— Vu,a«) of the particle 
with respect to the wall, and a stochastic force 



F' 



(16) 



which satisfies the fluctuation-dissipation relation and 
thus ensures the Boltzmann distribution to be the lo- 
cal equilibrium distribution. Here a is x, y, z and Xi,a is 
a Gaussian distributed random variable with mean zero 
and variance one: (x^.a) = 0, {xi,aXj,0) = SijSa/3- The 
viscous coupling between fluid and wall is realized by the 
locally varying viscosity j^uiLiz) with lolIz) = 1 — z/zc 
up to a cut-off distance Zc- Beyond z < Zc, the fluid 
satisfies the Stokes equation with an effective boundary 
condition of the form ([4]) The prefactor 7^ can be 

used to tune the strength of the friction force and hence 
the value of the slip length Sb- We note that is an ef- 
fective parameter and not related to the actual slip at the 
physical walls, which is usually nonzero even at Sb = 0. 
Within this approach it is possible to tune the hydrody- 
namic boundary condition systematically from full-slip 
to no-slip, and to derive an analytic expression for the 
slip length as a function of the model parameters [2^. 
Theorists tend to favor no-slip boundary conditions. In 
microchannels, howev er, p artial-slip boundary conditions 
are also observed 2^, 27, 2^. Here we will show results 
for both no-slip and partial-slip walls. 



Lattice-Boltzmann Approach 

In contrast to DPD, the Lattice-Boltzmann (LB) 
method can be seen as a discrete formulation of the 
Boltzmann equation on a lattice, which, by means of a 
Chapman-Enskog expansion leads to the Navier-Stokes 
equation in the incompressible limit [2^, [30| . The basic 
evolution equation for the mass density ni(r,t) assigned 
to a discrete velocity — Cia/r on a lattice node r at 
time t is given by 

nt{r+CiT,t+T) = ni{r,t) + ^Lij {nj{r,t) - nf{p,u)) , 

(17) 

where r is the time step, and is a set of vectors con- 
necting the grid point r to its nearest and next-nearest 
neighbors on a simple cubic lattice with lattice spacing a. 
The zero- velocity vector is also included (D3Q19 model). 
The collision matrix L relaxes the rii towards the local 
pseudo-equilibrium n^'^{p, u), which depends on the local 
mass density p{r,t) = X^i "»('*> ^) the local fluid ve- 
locity u(r, i) — (1/p) Xj ni(r, f)ci. The functional form 
of the pseudo-equilibrium is taken as the second-order 
approximation 



(p, u) = pa"' 1 



u • c, 



(u. 



2d 
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which maximizes the entropy of the underlying gener- 
alized lattice gas model [30|. Cg is the speed of sound 
and the parameters a^' are weight factors that depend 
on the magnitude of the velocity vectors but not their 
direction. For the D3Q19 model the respective values are 
flO = 1/3, fli = 1/18 and = 1/36, and the speed of 
sound is Cs = •\/l/3 a/r. In order to reproduce Brownian 
motion in a suspension, thermal fluctuations have to be 
introduced. The fluctuating lattice Boltzmann equation 
is given by 

ni{r + CiT,t + t) = n,(r,t) + 

{nj{Y,t)~-nf{p, u)) 



(19) 



For details regarding the implementation of the stochas- 
tic term n[{v,t) we refer to Ref. H^. A delicate task in 
lattice-based simulations is the correct coupling between 
the discrete nodes of the LB-solvent and the continuous 
positions of the (ion) particles. This can be achieved with 
a Stokes-likc friction force 1331 



■ 7nl 



= -C[V-u(R,t)]+f, 



(20) 



which is exerted on a solute particle moving through the 
surrounding viscous fluid. The random force f is required 
to guarantee thermal equilibrium of the coupled system, 
and its amplitude can be obtained from the fluctuation 
dissipation relation < fa{t)fp{t') >= 2(kBTSapS{t — t'). 
Momentum is exchanged between the particles and the 
LB fluid according to total momentum conservation in 
the system. The LB simulations were carried out for 
no-slip boundary conditions only, which were realized by 
applying bounce-back boundary conditions (s^ . 



Parameter Mapping 

Our scheme for mapping the parameters of the dif- 
ferent simulation methods onto each other is based on 
the requirement that the hydrodynamic flow phenomena 
should be the same. Therefore, the values of the parame- 
ters entering the Stokes equation ^ should be identical, 
i.e., the solvent density p and the dynamic viscosity rj. 
Matching the solvent density p is trivial, because it is an 
input parameter both in the DPD and the LB method. 
Matching the dynamic viscosity 77 is more difficult. This 
parameter is an input parameter in the LB method, but 
an a priori unknown intrinsic fluid property in DPD mod- 
els. A mean-field analysis of the DPD method [3^ shows 
that it mainly depends on the friction coefficient 7, the 
fluid density p, the temperature T, and the cut-off range 
of DPD interactions re- For given p and T, it can hence 
be adjusted by varying 7 and/or re- The dynamic vis- 
cosity can be determined either by measuring the auto- 



correlation function of the pressure tensor in free solu- 
tion, which is related to the viscosity via a Green-Kubo 
equation [3I] , or by simulating Plane Poiseuille flow in a 
confined microgeometry [2^ , respectively in free periodic 
boxes [s^l- The resulting flow profile is then fitted to a 
parabola Vx{z) = a{Zp — z^) with the magnitude 



2a 



(21) 



where F^: is the body force on the solvent which was 
applied to generate the Poiseuille flow. Having fixed p 
and 77, a third parameter remains which is related to the 
long-time dynamical behavior of single particles, i.e., to 
the self-diffusion. This is another intrinsic fluid property 
in DPD fluids. The LB method does not operate with 
particles originally, but in the hybrid LB/MD-mcthod 
(see section IILB), particles can be introduced which are 
coupled to the LB-nodes with a coupling constant C, (cf. 
Eq. (|20p ). Therefore the self-diffusion coefficient D can 
be tuned by varying the coupling constant ( until it is 
in agreement with the self-diffusion coefficient of a solute 
particle in the DPD simulations. This can be checked 
by comparing the velocity autocorrelation function of a 
particle in both systems which is connected to the self- 
diffusion coefficient 13 by a Green-Kubo expression, or by 
simple comparison of the mean-square displacement of a 
tracer particle [1^. 



SIMULATION DETAILS 

All simulations have been carried out with exten- 
sions of the freely available software package ESPResSo 
[45I . m, We use a cubic simulation box {12a x 

12tT X 12a) which is periodic in the x- and y— direction 
and confined by impermeable walls in the z direction. 
Electrostatics for homogeneously charged walls are cal- 
culated by P3M and the ELC-algorithm for 
2D + h slab geometries with an ELC gap size of 2a. 
See Refl4ll for a definition of the parameters. In ad- 
dition, the MMM2D algorithm \^ was used to 
study homogeneously charged walls in the same geom- 
etry. Ions have a mass m and carry a charge q. They in- 
teract via a Weeks-Chandler-Andersen(WCA)-potential 
44| Vij = 4e{{a/r)^'^ — {a/r)^) + e with energy parame- 
ter e and cut-off distance rwcA = 2^/^cr, and a Coulomb 
potential with Bjerrum length Xb ~ 1.0a. The param- 
eters e, a, and m shall be used as the natural units of 
our system hereafter. To study the counterion-induced 
electroosmotic flow, the walls were charged by placing 
discrete charges q randomly all over them. The charge 
conflguration in each wall was identical for both methods 
with a surface ion density as (the surface charge density 
is then given by a a = qcs)- The number of counterions 
in the fluid matches the number of charges in the wall. 
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TABLE I: Charge q, surface ion density as, Gouy-Chapman 
length ^ and coupling constant H for different electrostatic 
coupling regimes. 



Regime 










Weak coupling 


1 


0.2083 


0.764 


1.31 


Intermediate coupling 


2 


0.042 


0.955 


4.19 



such that the system is overall neutral in the case of inho- 
mogeneously charged walls. The case of homogeneously 
charged walls can be realized simply by having both walls 
uncharged. This is due to the fact that the electric fields 
generated by the two plates cancel each other exactly 
within the slit, and the charge density profile is gener- 
ated only by the internal Coulomb interaction between 
counterions. 

Different coupling regimes were obtained by varying 
the ion charge g, the surface ion density CTs and the cor- 
responding counterion number from 12 to 60. The pa- 
rameters are shown in Table |I] together with the resulting 
coupling constant 2. The density of solvent particles was 
chosen pi = 5.75a~^, and the temperature k^T = le in 
all simulations. 



Dissipative particle dynamics simulations 

In the DPD simulations, the walls are located at 
Zwaii = (0, 10)(T. They interact with fluid particles and 



ions via a WCA-potential [4J] with the same parameters 
as above {ewaii = ^,<ywau = cr). The solvent particles 
have the same mass as the ions (m), but do not interact 
with other particles except the walls. The DPD fric- 
tion coefficient was chosen 7 = 5.0(T^^(TOe)-'^/^ and the 
cutoff range of the DPD interactions Tc = l.Ocr. Tun- 
able slip boundary conditions [2^ were used with fric- 
tion coefficients 7^ = 0.96, 1.4049, 3.1(7-Hme)i/2 and 
7i = 6.1cr~^(me)^/^ with the friction range Zc = 2.0(7 



starting at z^uaii = 
5t = 0.01cr(m/e)i/2. 



(0,10)cr. The DPD timestep is 



Lattice Boltzmann simulations 

The LB simulations were carried out using the D3Q19 
model with 24^ solvent nodes. The walls for the ions are 
placed as in the DPD simulations at z^aii = (0, 10)cr. 
Since the zero-velocity surface in a LB simulation is lo- 
cated, at first order in viscosity, halfway in between two 
rows of nodes, a shift has been added to the WCA wall 
potential in order to match the position of the zero- 
potential and zero-velocity planes. The grid spacing of 
the LB fluid is a = O.Scr. The coupling constant of the 
fluid with the ions is C = 1.98tT-i(me)V2 which was de- 
rived by the mapping scheme proposed in Section 5. The 
integration timestep for the ions as well as for the fluid 



TABLE IL Time needed for computing a single integration 
step in the DPD method with tunable-slip boundary con- 
ditions (DPD-hTSC) and electrostatics (DPD+TSC-(-CS) in 
comparison to the LB method with bounce-back boundary 
conditions (LB-I-BBC) and electrostatics (LB-I-BBC-I-CS). 



Methods 


DPD 


LB 


DPD-hCS 


LB-I-CS 


Time/step [s] 


0.09 


0.01 


0.22 


0.14 



was 5t = T = 0.01(T(m/e)"'^/^, bounce-back boundary con- 
ditions were applied on the fluid at the wall positions to 
create no-slip boundary conditions. 



RESULTS 
Computational cost 

An important criterion when choosing a mesoscopic 
simulation method is its computational cost and effi- 
ciency. We have measured the time to calculate a sin- 
gle time step in both methods on an Athlon© MP2200+ 
CPU. The values are presented in Table [ill The first 
two columns show the time spent on an uncharged sys- 
tem with 4320 solvent particles (DPD) with tunable-slip 
boundary conditions (TSC) or 1728 solvent nodes (LB) 
in the above mentioned microgeometry with bounce- 
back boundary conditions (BBC). The LB-method is nine 
times faster than the DPD-method. We should note, 
however, that these values strongly depend on the density 
of solvent particles, or the number of solvent nodes, re- 
spectively. The last two columns show the corresponding 
values for a charged system (CS) with the P3M in com- 
bination with the ELC algorithm for discretely charged 
walls. It is evident that most of the time is spent on the 
calculation of the electrostatic interactions for 60 ions 
and 60 counterions [as = 0.208(t~^). If electrostatic in- 
teractions are considered, both methods are comparable 
with respect to their computational cost and efficiency, 
although it has to be noticed that the efficiency of P3M 
strongly depends on the chosen parameter values (Ewald 
parameter a = 2.1875, mesh size 32'^). Depending on the 
choice of parameter and of system size the computational 
load of computing electrostatic calculation can be even 
lower than that of hydrodynamics. 



Fluid properties 

We first consider the dynamic properties of the bulk 
fiuid. To measure the dynamic viscosity of the DPD fluid, 
which is needed for the parameter mapping, we fitted 
a plane Poiseuille flow as in Ref. [H. This procedure 
yields the value rj = (1.334 ± 0.003)(T-2(me)i/2 both for 
uncharged and charged fluids. Next we calculated the 
effective diffusion coefficient D for a single tracer particle 
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FIG. 1: Normalized velocity autocorrelation function for 
a DPD-fluid particle (open circles) with number density 
p = 3.75cr~^ DPD friction coefficient 7 = 5.0cr"^(me)^''^ 
and for a LB/MD-fluid particle (filled triangles) with the 
same density and coupling constant = 1.98a~^ {meY^^ 
The characteristic decay time for the DPD-method is 
TDPD = (0.5162 ± 0.0008)(T(m/e)^/2 and for the LB-method 
TLB = (0.5218 ± 0.0006)(T(m/e)^/^ Inset: Mean square dis- 
placement for a fluid DPD particle (open circles) and for a 
coupled LB/MD particle (flUed triangles) compared to Eg. 1241 
(solid line). 



in the DPD-mcthod and matched it to the tracer diffusion 
in the LB-modcl. The self-diffusion coefficient D can 
be obtained by exploiting a Grcen-Kubo expression [38| . 
which relates D to the velocity autocorrelation function 
(see Fig. H]) 



D=l- f dt< v(t)v(io) > 

J Jto 



(22) 



Fig. [T] shows that the velocity autocorrelation function 
decays exponentially over at least two orders of mag- 
nitude. Approximating the integrand in Eq. (|22p by 
this exponential law we obtain D_dp_d = (0.2581 ± 
O.GOG4)cr(m/e)-i/2 for the particle in the DPD fluid (at 
late times, one expects to see an algebraic long-time tail 
in the velocity autocorrelation function, but its contri- 
bution to D is small and well within the error). A 
coinciding value of the self-diffusion coefficient in the 
LB/MD-method can be obtained by setting the coupling 
constant to C = 1.98cr~^(TO/e)^/^. The two correspond- 
ing velocity autocorrelation function then closely match 
each other (Fig. [T]) and the resulting self-diffusion co- 
efficient for LB/MD-fluid particles is given by Dlb = 
(0.2609 ± 0.0003)cr(m/e)-i/2. 

Alternatively, the self-diffusion coefficient can be de- 
termined directly from the mean-square displacement of 
a single solvent particle at late times 



D = hm 

t—*oo 



< (r,(t)-r,(^o))^ > 
6t 



TABLE III: Slip lengths for different layer friction coefficients 
7l, determined for uncharged fluid flow. The numerical re- 
sults are compared to the theoretical prediction of Ref.25. 



7i 


0.96 


1.4049 


3.1 


6.1 


cm 


1.399 ±0.385 


0.782 ± 0.246 


0.248 ±0.231 


0.000 ±0.197 




1.249 


0.780 


0.226 


0.000 



The results for the mean-square displacement are shown 
in the inset of Fig. [TJ In agreement with standard the- 
ories [3 IS], ballistic behavior (~ t'^) is observed at 
short times {t < 0.75CT(m/e)^/^), but diffusive behavior 
dominates after a characteristic time which is roughly 
t > 10cr(m/e)i/2 for our model parameters. A lin- 
ear regression of the diffusive regime yields Dupjj = 
(0.2698 ± 0.0002)cr(?7i/e)-i/2 ^nd Dlb = (0.2617 ± 
0.0005)cr(m/e)~^/^. which is in rough agreement with the 
Green-Kubo results. It is remarkable that the numerical 
results obtained with both methods lie almost on top of 
each other. Assuming that the velocity of a single particle 
propagates according to an Ornstein-Uhlenbeck process, 
the full-time mean-square displacement can be calculated 
analytically as a function of the effective friction coeffi- 
cient = ksT/D and is given by 



< (r,(t)-r,(to))' > 



6 — —t H — (1 



mQ 



(3- 



). (24) 



(23) 



This prediction is shown as the straight line in the inset of 
Fig. [TJ It is in reasonable agreement with the numerical 
results. 

Finally, we consider uncharged DPD fluids in slit ge- 
ometry and establish the values of the slip length 5b and 
the hydrodynamic boundary positions zb for different 
values of the surface friction parameter 7/,. By a com- 
bination of plane Poiseuille flow and plane Couette flow 
simulations it is possible to determine 5b and zb indepen- 
dently [25I . The resulting hydrodynamic boundary posi- 
tions were found to be located at \zb\ = (3.866 ±0.265)(7, 
and the corresponding results of the slip length compared 
to the theoretical prediction of the analytical expression 
derived in [2^ are shown in Table IIIIl 

Having matched the model parameters in the DPD and 
the LB/MD method, we can now proceed to simulate the 
counterion-induced electroosmotic flow and to compare 
the results. 



EOF: Weak-coupling regime 

First we consider the counterion distribution in the 
channel with inhomogeneously charged walls in the weak- 
coupling regime (S = 1.31). DPD simulations were car- 
ried out for ten different external electric field strengths 
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FIG. 2: Counterion density distribution for the DPD-fluid 
in the weak-coupling limit (H = 1.31) between two charged 
walls for external electric field strengths between Ex = 0.1 — 
l-QkaT / ea. The external perpendicular electric field does 
not perturb the ion density distribution. The dashed line 
presents the theoretical prediction of the Poisson-Boltzmann- 
Theory for an ion density in the middle of the channel of 
po = (0.0176 ± 0.0001)(J~''. Inset: Comparison of counterion 
distribution for the DPD-method (circles) and for the LB- 
method (triangles) with a field strength = l.OfcsT/ecr. 
The straight line is again the theoretical prediction of the 
Poisson-Boltzmann- Theory. 



between = 0.1 - l.Ofcsr/ecr. Fig. H shows that 
the counterion distribution between the walls is not 
perturbed by the applied fields nor the discreteness of 
charges in the walls in the accessible z-iaxige. For all 
fields, it is in very good agreement with the theoretical 
prediction of the Poisson-Boltzmann equation ([7]) with 
the only fit parameter po = (0.0176 ± 0.0001)cr"3^ The 
inset of Fig.[2]compares the counterion distribution in the 
DPD-fluid with the corresponding ion distribution in the 
LB/MD-fluid for a field strength = IMsT/ea. The 
simulation results are nearly identical and follow again 
closely the prediction of the electrostatic theory. 

In addition, we also measured the electric field E{z) 
by a test charge method. The results are shown in 
Fig. [21 they also agree very well with the prediction of 
the Poisson-Boltzmann theory, Eq. ([8]) . 
In addition, we also measured the electric held E{z) by a 
test charge method. The results are shown in Fig. [31 they 
also agree very well with the prediction of the Poisson- 
Boltzmann theory. Eq. ([5]). 

The influence of the partial-slip boundary conditions 
is presented in Fig. [H for varying field strength E^- In 
agreement with the theoretical prediction, Eq. ([9]), we 
find that oc E^ for all slip lengths. Varying the slip 
length has a quite dramatic effect on the magnitude of 
the fiow profiles. All numerical results are in reasonable 
agreement with Eq. especially if one bears in mind 
that the theoretical curves have an uncertainty due to 
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FIG. 3: Electric field E{z) in the solvent measured by a test 
charge method in the weak-coupling limit (H = 1.31). The 
dashed line shows the theoretical prediction of Eq. (|8]). 




FIG. 4: Flow profiles for the DPD- method with various field 
strengths Ex = 0.8 — l.OfcsT/eo- for varying slip lengths 
in the weak coupling regime (H = 1.31). The hydrody- 
namic boundary positions for the DPD-method are at \zb\ = 
(3.866 ± 0.265)a. The straight lines represent the theoretical 
prediction of Eqn.([9ll. 



the numerical error of 5b- 

The comparison of the DPD- and the LB/MD flow 
profiles is finally presented in Fig. [5] for no-slip boundary 
conditions. Here we have shifted the effective channel 
width in the DPD method from 8cr to 8.28(7, such that 
the hydrodynamic boundaries were now located at \zb\ = 
4.030±0.357cr like in the LB/MD-fluid. Fig.[5lshows that 
the flow profiles obtained with both methods are then 
identical and agree well with the theoretical prediction, 
Eq. evaluated with 6b = and \zb\ = 4(T. 
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0.15 




0.05 



FIG. 5: Flow profile for the DPD-(circles) and the LB/MD- 
method (triangles) with no-slip boundary conditions. The 
straight line is the theoretical prediction of Eqn. @ with 
\zb \ = 4.0a and 5b = O.Ocr. 




FIG. 6: Flow profile for various slip lengths with homoge- 
neously (filled diamonds) and inhomogeneously charged walls 
(circles) for Ex ~ LOftsT/ecr in comparison to the theoret- 
ical prediction of Eqn. @. Inset: Counterion distribution 
with (TA = 0.208e(T~^ between homogeneously respectively 
inhomogeneously charged walls. Only slight deviations can 
be observed. 



EOF: Homogeneously and Inhomogeneously charged 

walls 



The results shown up to now have been obtained in 
system where discrete charges were placed randomly in 
the walls. For comparison, we have also studied sys- 
tems with homogeneously charged walls with DPD simu- 
lations. Surface charge inhomogeneities may lead to elec- 
trofriction and slow down the fluid at the wall [sot . Since 
the electrostatic potential between two equally and ho- 
mogeneously charged walls is constant, the charges on the 
wall can be omitted and it suffices to study a fluid with 
ions confined between uncharged walls. The simulations 



were carried out using the MMM2D-algorithm [43, |43[ 
and the parameters given in Table IIIII for the weak- 
coupling (Poisson-Boltzmann) regime. The ion density p 
was 0.0525cr~'^ in all simulations. Fig. [5] presents the nu- 
merical results for various slip lengths for homogeneously 
and inhomogeneously charged walls, compared to the 
theory of Eq. ([9]) for an external field of = LOfesT/ecr. 
All flow profiles for homogeneously respectively inhomo- 
geneously charged walls are identical. These results are 
consistent with theories [H^, [5l| - the electrostatic inter- 
action is moderate and drastic deviations have only been 
reported for strongly interacting systems [H^l- In the 
weak-coupling regime, the infiuence of electrofriction on 
the flow profiles both for no-slip and partial-slip bound- 
aries is thus negligible. 



EOF: Intermediate-coupling regime 

Finally, we consider a situation where no analytical 
theory is available, and study the intermediate coupling 
regime with S = 4.19. The simulation results for the 
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FIG. 7; Counterion density distribution with p — 0.0104(t~ 
in the intermediate coupling regime with H = 4.19 for a sur- 
face ion density a a — 0.042(7^^. The straight line shows the 
theoretical prediction of the strong-coupling limit (Ea. (|10|) '). 
Neither the strong coupling theory (straight line) nor the 
Poisson-Boltzmann theory (dashed line at the bottom) are 
able to reproduce the results. The fit function is represented 
by the dotted line. Inset: Blowup of the z-range between 
—4.0 and —2. Oct in a logarithmical plot. The counterion 
density distribution decays exponentially proportional to the 
Gouy- Chapman length of /i = 0.955cr. 

counterion density distribution are presented in Fig. [T] 
The distribution can neither be fully described by the 
strong coupling theory (Eqn. (fTO|) ). shown as straight 
line, nor by the Poisson-Boltzmann-Theory (Eqn. (O), 
shown as dashed line. The inset of Fig. [7] shows that 
the counterion distribution close to the walls decays ex- 
ponentially with the Gouy-Chapman length /i, indicat- 
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FIG. 8: Counterion induced electroosmotic flow for the in- 
termediate regime (H = 4.19) with p = 0.0104(j~^, E^; = 
lO-OfesT/eo- and surface ion density as = 0.042(7~^ with 
7i = 6.1cr-i(me)^/2 {5b = (0.000 ± 0.197)a) and 7i = 
6.1cr-^(me)^''2 (Sb = (0.248 ± 0.231)o-)for homogeneously 
(filled diamonds) and inhomogeneously charged walls (cir- 
cles). The straight lines show the calculated flow profile for 
|2:s| = 3.866cr and for slip lengths as indicated. 



ing the presence of a layer of highly adsorbed counte- 
rions. This corresponds to the behavior predicted by 
the strong-coupling-Theory (Eqn. (fTUl) ). In the bulk, 
however, the data are better described by the Poisson- 
Boltzmann theory in agreement with Ref. [2^ . Thus, 
the intermediate regime bears characteristics of both the 
strong-coupling and the Poisson-Boltzmann regime. For 
the former regime, a recently developed method (s^ 
might turn out to be useful. These results are consistent 
with Refs. 22, 5^, where it was shown that in case of a sin- 
gle charged plate, neither the Poisson-Boltzmann nor the 
strong-coupling-theory are applicable at distances from 
the plate d' in the range Ve < d'/M < Inserting the 
parameters of S and fj,, this corresponds to the z-range 

\z\ > 2(7. 

Although we have no analytical theory, we can still fit 
the counterion density with a purely heuristic function 
whose functional form is inspired by the prediction of 
the strong coupling theory, Eq. (jlOp . 



p\z) ^ (e-^^-^n)/, + + /corr(z), (25) 

where p*-^^ is a fitting parameter and /corr(^) a phe- 
nomenological correction function. A good fit can be 
obtained with /corrC^) = P*2) cos((/)z) with and p*2) 
being two additional fit parameters. The result of the fit 
to this expression (with fitting parameters values — 
(8.31 ± 0.10) • 10-6(7-3, = (3-18 ± 0.05) • lO-^cr-s, 
and (j) = (0.416 ± 0.022)cr--'^, respectively) is presented 
as a dotted line in Fig. [T] Combining Eq. (pS)) with the 
Stokes equation, Eq. and integrating it using Eq.Q 
as boundary condition, allows to reproduce the EOF pro- 



files very nicely (Fig. [8]). Our result demonstrates that 
the Stokes equation is still applicable in the intermediate 
coupling regime. Furthermore, Fig. [5] shows that we still 
do not observe a noticeable effect of the surface charge 
inhomogeneity, compared to simulations with homoge- 
neously charged walls. Thus, electrofriction plays a neg- 
ligible role in the intermediate as well as in the weak- 
coupling regime and seems to be a prominent effect only 
in the strong-coupling regime. 



SUMMARY 

We have carried out two types of mesoscopic sim- 
ulations, namely, Dissipative Particle Dynamics and 
coupled Lattice-Boltzmann/Molecular Dynamics, of the 
countcrion-induced EOF in slit channels. We have pro- 
posed a mapping scheme that allows to match quanti- 
tatively the parameters of the two models. The cost of 
calculation time for electrostatic problems in both meth- 
ods turned out to be nearly identical. 

We have considered different electrostatic coupling 
regimes. In the weak coupling regime, the numerical re- 
sults for the ion density and the electric field are in good 
agreement with the predictions of the Poisson-Boltzmann 
theory. We have derived analytical expressions for the 
corresponding EOF profile in presence of no-slip as well 
as partial-slip boundary conditions, which also agree well 
with the simulations. In addition, we have considered the 
intermediate coupling regime where no analytical theory 
is available. A heuristic function was used to fit the coun- 
terion distribution in this transient regime. The Stokes 
equation can still be used to calculate the corresponding 
EOF profiles, with very good results as compared to the 
simulations. 

In all systems under consideration, electrofriction ef- 
fects were found to be negligible. In contrast, the pres- 
ence of partial slip changes the magnitude of the flow 
profiles drastically even if the slip lengths are very small 
[3] . This may facilitate the generation of flow profiles in 
microfluidic applications. 
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